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Abstract The increasing number and variety of extrasolar planets illustrates 
the importance of characterizing planetary perturbations. Planetary orbits are 
typically described by physically intuitive orbital elements. Here, we explic- 

■ itly express the equations of motion of the unaveraged perturbed two-body 
problem in terms of planetary orbital elements by using a generalized form of 
Gauss' equations. We consider a varied set of position and velocity-dependent 

P I ■ perturbations, and also derive relevant specific cases of the equations: when 

. they are averaged over fast variables (the "adiabatic" approximation) , and in 

(-H I the prograde and retrograde planar cases. In each instance, we delineate the 

^ ^ properties of the equations. As brief demonstrations of potential applications, 

we consider the effect of Galactic tides. We measure the effect on the widest- 
known exoplanet orbit, Sedna-like objects, and distant scattered disk objects, 

■ particularly with regard to where the adiabatic approximation breaks down. 
. The Mathematica code which can help derive the equations of motion for a 

user-defined perturbation is freely available upon request. 



Keywords Perturbation Methods; Computer Methods; Planetary Systems; 
Comets and Meteors 
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\^ • 1 Introduction 

■ The discovery of extrasolar planets challenged our understanding of the forma- 

I tion and subsequent evolution of planetary systems. The variety of dynamical 
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architectures den ionstrated by pulsar pl anets (e.g. Wolszczan fc Frail. 19921) 



Hot Jupitcrs ( e.g. Mayor fc QuelozL 1995 ). wide-orbit pl anets (e.g 



Zapatero Osorio et al.l.l200C ; Sumi et all. 1201 ll). highly inclined and retrograde 



Kuzuhara et al 



2011; Luhma n et all. 2011 ) free-floating plan ets (e.g. Lucas fc RocheL 2000l : 



planets fe.g. lAnderson et al.l . l2010l:lTriaud et aLlbOldlWinn et al.ll20ldHKaib et al 
I2OII), closely-packed planets (e.g. Lissauer et al.f 2011 ). circumbinary planets 
(e.g. Sigurdsson et al.L 2003 : Dovle et al.l . 201ir and the Solar System defy a 
single, simple explanation for their formation. 

After formation, these systems are subjected to internal forces, from plan- 
ets and smaller bodies, and external forces, such as the singular or repeated 
local close encounters with individual passing stars, the Galactic tide, tidal 
tails, molecular clouds, globular clusters or even dark matter substructures. 
Further, a star's passage into and out of spiral arms, and the deformation of 
the Galactic tide due to close encounters or collisions with other galaxies can 
affect orbiting planetfl One known exoplanet is thought to be of extragalactic 



origin , perhaps originating from a disrupted satellite galaxy ([Setiawan et al 



2OIOI ) . The vast populat i on of free-floating planets - nearly two per main se- 
quence star ( Suini et"all 201 1|) - which ca nnot be explained by planet-planet 
scattering alone (jVeras fc RavmondLlioli ). perhaps demonstrates that exter- 
nal disruption of planetary systems is an endemic feature of the Galactic disk. 

Some of these external forces may be modeled as perturbations on a two- 
body star-planet system. Although the two-body problem is solvable analyt- 
ically, the perturbed two-body problem is generally not so. Nevertheless, ex- 
pressing these equations of motion entirely in terms of planetary orbital ele- 
ments garners intuition for how planets are affected by external perturbations 
and helps one obtain analytical solutions in limiting cases. Here, we present a 
method for performing this conversion for arbitrary perturbations. 

Suppose the equations of motion for a single planet orbiting a single star 
are subject to extra accelerations T. Then the equations of motion can be 
expressed as: 
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■T^.x + T^yy + Tuu^+Tuy— (1) 



dz 

' '^zz^ ~l~ '^ww~J^ (3) 

where m-^ and nip represent the masses of the star and planet, respectively. 
In the context of a single planet orbiting the Galactic center, we may take 
{x,y,z) to be a non-rotating rectangular coordinate system centered on the 



^ The Milky Way and Andromeda will collide i n ~ 2-5 Gyr, before many Galactic stars , 
inclu ding the Sun, turn off of the main-sequence lICox &: Loebl . |2008| : Ivan der Marel et al.L 
120121) . 
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star, as in iHeisler fc Tremainel (1986). The primary goal of this paper is to 
express Eqs. ([I])-© analytically in terms of the time evolution of planetary 
orbital elements. The secondary goal is to present the method by which one 
can repeat the procedure for other perturbational forms not included in Eqs. 

m-m- 

In Section 2, we provide background for the perturbed two-body problem 
and present the algorithm used to derive the equations. We then present the 
general equations in Section 3. Sections 4 and 5 specialise to the adiabatic 
and planar adiabatic cases, for which substantial simplifications are possible. 
We provide applications to perturbations derived from the Galactic tides in 
Section 6 and summarize our results in Section 7. 



2 Background and Derivation Algorithm 

In order to describe fully the position and velocity of a planet, 6 orbital ele- 
ments must be employed. Four often-used elements are the semimajor axis, a, 
eccentricity, e, inclination, i, and the longitude of ascending node, ]7. A fifth 
orbital element usually includes the longitude of pericenter, zu, or the argu- 
ment of pericenter, lu. The difference is important here, and they are related 
through oj + n = w. All five of these elements fully describe the shape of a 
planet moving on a Keplerian ellipse. In the absence of any other planets or 
forces, these elements remain fixed, as in the classic two-body problem. The 
sixth element, which provides the location of the planet along its orbit, can 
be described by one of several different measures, including the mean anomaly 
M, true anomaly /, eccentric anomaly E, true longitude 9, mean longitude 
A, or argument of latitude v. We adopt / as our sixth element due to its in- 
tuitive geometric interpretation and with an eye towards future studies that 
might wish to focus on planets residing at locations close to either pericenter 
or apocenter. 

The perturbed two-body problem has been a subject of extensive study 
and is applicable to several fields of astrophysics. The work of iBurn j (1976), 



subsequently popularised bv lMurrav fc Dermotd (jl999l pp. 54-57), provides a 
mechanism to obtain analytic equations for da/dt, de/dt, di/dt, and dQ/dt 
arising from a small perturbative force with a given prescription for radial, 
tangential and normal components. This line of attack can be traced back 
ultimately to Gauss (see e.g., section 9.13 of Brouwer & Clemence 1961). 

However, there is an alternative. Recently, a re-analysis of the derivation of 
Lagrange's Planetary Equations, which describe the perturbation from a third 
body, revealed that a generalized gauge may be adopted in the derivation. 
The consequences and extensive description of this "generalized gauge theory" 
is provided in Efroimskv & Goldrcich (2003, 2004), Gurfil (2004), Efroimsk^ 
(t2005allbl . l2006l) . lGurfill (l2007t) and lGurfil fc BelvaninI (|2008l) . This theory also 
yields analytic equations for da/dt, de/dt, di/dt, and dQ/dt but also directly 
for dM/dt, arising from a perturbative force with a given prescription for 
each Cartesian component. In principle, either approach can be used for the 
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purposes of this paper; we use the latter to derive the equations of motion due 
to its generality and compactness. 

To derive the equations of motion in orbital elements for any perturbing 
acceler ation, we choose the fo l lowing form, similar to Eq. (22) of lEfroimskv 
(|2nn5ah and Eq. (16) of lGurfill (|2n07h : 



d£ 

dt 



= Mi 



(4) 



where /3 = {a, e, i, Q, uj, A/q}, the subscript "0" refers to the initial value, and 
* is the gauge. 

The AA represent matrices that are expressed entirely in terms of orbital 
elements, and Z\A is the acceleration caused by a perturbative force on the 
system. Although Z\A is typically written and referred to as a force, the quan- 
tity actually represents an acceleration. Equation Q is useful because all three 
matrices are system-independent, and can be prccomputed before any pertur- 
bative force or gauge is applied. 

Assume A is expressed as a column vector. The matrix A^i is then the 
transpose of the "Poisson Matrix", which is composed of Poisson Brackets, 
and is the negative inverse of the "Lagrange Matrix", which is composed of 
the Lagrange Brackets. This matrix reads: 



Ml 









2 

na 







VT- 



- \ 

n.n. 



na^ sin 

1 



nea^ 
1-e 



na^ tan i^l — e^ 





na^ sin 







na^ tan i\/l~e^ 







where the orbital elements associated with each row starting at the top are a, 
e, i, fi, Lo and Mq, respectively. Note that these matrix entries represent the 
coefficients in the appropriate form of Lagrange's Planetary Equations. The 
matrices M.2 and M.3 are composed of partial derivatives of r and v, namely: 
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In order to compute these partial derivatives, we consider the relation between 
r and v in a fiducial inertial reference fr ame. These relations can be found 
in classical mechanical textbooks such as iTaS (|l985l pp. 35-36). A standard 
choice is: 
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7^ 



/ Q^l-£)_COS/ \ 

1+e cos / 
a^l — e^^ sin / 
l+e cos / 











= n 







-na sin / 




(5) 



where TZ is the matrix 

(cos I? cos w — sin ]7 sin cos i — cos^?sina; — sinJ7cosa;cosi sin i? sin i 
sin ]7 cos + cos J7 sin a; cos i — sini7sinw + cosJ7cosa;cosi — cos^?sini 
sin uj sin i cos oj sin i cos i 

(6) 

Our choice of fiducial reference frame sets the x-axis to be along the major axis 
of the ellipse, where the pericenter is in the positive direction. The matrices 
TMi, A^2, -A^s are independent of the system being studied. Their computation 
does, however, contain some subtleties which are worth recording: i) the mean 
motion n in 3 is a function of a, ii) the true anomaly / is a function of 
the time t as well as Mq, e and a, iii) by definition, Mq is defined at i = to- 
Therefore, we need to compute the partial derivatives with respect to the true 
anomaly: 



df df BE 



dMo dE dMo 



l + ecos/\ /l + ecos/\ (1 + e cos/) 



(1 



^3/2 



(7) 



where the expression for the first partial derivative is fromlBroucki (|l970l ) - 
who also provides tables of partial derivatives, Lagrange Brackets and Poisson 
Brackets for the unperturbed two-body problem - and the second is from 
Kepler's equation and from the relations between eccentric anomaly E and 
the true anomaly / H: 



cos£^ — 



cos/ 



1 + e cos / ' 



sini? = 



sin /\/l — e^ 
1 + e cos / 



(8) 



Further, we have: 
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(2 + ecos/) (9) 



and 



df _ df dE _ f l + ecosf \ f 3n{t~to) \ _ 3 ^, (l + ecos/)' 

da- dEda^[ J\ 2r J 2^^' "'ail-e^f 

(10) 



2 Note that Eq. (Ill) of lEfroimsk"^ ll2005ar) for dM/df contains a typo. 
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where the partial derivatives ol E with respect to e and a are from iBroucke 
(|l970l ). Importantly, one cannot impose the condition t = to until after M2 
and M3 have been computed. 

The time evolution with respect to E 01 f is often more desirable to in- 
vestigators than the time evolution of Mq, as is true in this study. We can 
convert one equation into another by taking the total time derivative of Ke- 
pler's equation evaluated a.t t ~ to, where M — n{t — to) + Mq. The result 
is: 



dE 



1 



1 — e cos E 



dMo 
dt 



smE — 
dt 



(11) 



We note that the terms with the time derivatives of the semimajor axis vanish. 
Then wc finally obtain: 



dt 



1 
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fl 



e cos 
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sin/ 



(12) 



We choose $ = in order to obtain a simply interpreted set of resulting 
orbital elements. However, there may be instances when choosing a non-zero 
gauge provides an insight or simplification to the problem considered. 

As a check on these equations, we used them to re-derive t he equations 



of mo tion for the two-body problem perturbed by mass loss in IVeras et al 
(120 Uh . 



3 General Equations 

The equations below are well-suited for describing a planet on an arbitrarily 
wide and inclined orbit. 



3.1 Orbital Element Equations 
We apply: 



AA = I TyxX + TyyU + T^^i -I- T^^y 



(13) 



to the formalism in Section 2 in order to derive the equations of motion in 
terms of orbital elements: 
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— {TyxCs — Tyyd) {Ci COS n cos i — C2 sin Q} 
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2a 



+ ^uu {Ci sin n cos i + C2 cos Q}^ + Tyy {Ci cos i7 cos i ~ C2 sin H}^ 
— O^uv + "^vu) {Ci sin f? cos i + C2 cos i7} {Ci cos cos i — C2 sin J7} 



(14) 



(l-e^)» 



at 2n(l + ecos/) L 
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2 (1 + ecos/j [ 

+ {Ci sin J7 cos i + C2 cos f2} {Ce sin J7 cos i + C5 cos 12} 

+ Tyy {Ci cos i7 cos i — C2 sin i?} {Ce cos Q cos i — C5 sin f2} 

— Tuy {Ci COS f2 cos i — C2 sin i?} {Ce sin i? cos i + C5 cos /?} 

+ T„„ {Ci sin J? cos i + C2 cos f2} {—Cq cos i7 cos i + C5 sin 



(15) 



(1 — e^) ^ sini 



Tzz {cos i cos (/ + oj) sin (/ + w)} 



dt ~ n{l + ecosff 

- {'^xxCs - T^^yCi) {sin Q cos (/ + w)} + (T^a^Cs - Tyyd) {cos cos (/ + w)} 
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H 7 — i^^ {Ci cos?} 

1 + e cos / [ 

— T„„ sin i7 {Ci sin J7 cos i + C2 cos i7} + Tyy cos ^? {— Ci cos fl cos i + C2 sin i7} 

+ Tuv sin J7 {Ci cos fl cos i — C2 sin Q} + cos J7 {Ci sin i7 cos i + C2 cos Q} (16) 
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— Tuu sin Q {C\ sin i? cos i + C2 cos /?} + cos i? {— Ci cos i? cos i + C2 sin J?} 

+ T„„ sin Q {Ci cos i? cos i — C2 sin J?} + cos i? {Ci sin f2 cos i + C2 cos J?} (17) 



dui 
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e sin (/ + w) + -Cg sin^ i 
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— iXxxC^ — TxyCi) {Cs sin J7 cos i — C7 cos i7} 
+ {TyxCs — Tyyd) {Cs cos n cos i + C-j sin i7} 

+ TT—r, — ^ ?r Ci ^t«M) I — 2e (cos / sin w + cos cj sin / cos^ 

2e (1 + ecos/j L 

+ sin^ i (sin (2/ + cj) — 3 sincj) } 

+ 3^M« {Ci sin i7 cos i + C2 cos fl\ {— Cs sin i7 cos i + C-j cos J7} 

— Ti,„ {Ci cos J7 cos i — C2 sin 12} {Cs cos J7 cos i + C7 sin ^2} 

— "fuv {Ci cos J7 cos i — C2 sin ^2} {— Cs sin fl cos i + C7 cos SI} 

+ '^vu {Ci sin n cos 1 + 6*2 cos i7} {Cs cos J7 cos z + C7 sin SI} 



(18) 



df n(l + ecos/) 

~ (l-e2)^/^ 2en(l + ecos/) 



{sin^ i sin (/ + cj) [Cg + 2e sin (/ + w)] } 
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^ / 7T ^u^u. {C1C9 sin^ i} 

2e (1 + ecos/j [ ^ 

3^jju {Ci sin fl cos i + C2 cos ^2} {— Cg sin fl cos i + C7 cos Q} 
Tyy {Ci cos J7 COS i — C2 sin 12} {Cg cos SI cos i + C7 sin S2} 
T^uv {Ci cos 12 cos i — C2 sin i7} {— Cg sin cos i + C7 cos SI} 

Tvu {Ci sin /2 cos i + C2 cos S2} {Cg cos SI cos « + C7 sin S2} 



dio .dS7 rt(l + ecos/) 



(l-e2 

where the quantities Ci , . . . Cg are given in Appendix A. Ah C terms are on 
the order of unity, and so the only variables in the square brackets not of the 
order of unity are the T values. 

Note the remarkably simple form taken on by df /dt in Eq. ([20)) . where the 
last term represents the two-body term. This form allows us to quickly derive 
some other results. The argument of latitude w = / + cj is a quantity which 
evolves with time according to: 



(19) 
(20) 



dv 
It 



. dS2 n (1 + e cos /) 



dt 



(l-e2) 
u + ^2 is 



3/2 



The evolution of the true longitude 

dO , . dSl n(l + ecosf)^ 



(21) 



(22) 
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3.2 Properties 

Equations (|14p - are completely equivalent to Eqs. (H} - (12I)- However, 
Eqs. - let us deduce properties of planetary motion, which may not 
otherwise be apparent: 

First, note that all the T terms are decoupled, allowing us to easily consider 
special-case physical situations by striking out the relevant terms. The rates of 
change of eccentricity, inclination, longitude of ascending node and longitude 
of pericenter (|(ie/(it|, \dflldt\ and are all oc 1/n cx c?^"^ for the 

position-only T terms, whereas the rate of change of semimajor axis \dal dt\ cx 
ajn (X a^/^. Turning to velocity- based T terms, |(ie/di|, |di/d<|, \dflldt\ and 
\dioldt\ are all independent of a, whereas \daldt\ oc a. 

We can also deduce some properties of the evolution of special orbits. 
Planar orbits (i = 0) will remain planar, but circular orbits will not remain 
circular. Likewise, polar orbits will not remain polar. Further, the sign of dijdt 
will be different depending on whether i = 90° or i = 270°. 

As \duj /dt\ oc 1/e and \df /dt\ oc 1/e, then for circular orbits, neither ui nor / 
are well-defined physically nor mathematically. Although Q is not physically 
defined for planar orbits, fi still remains well-defined mathematically, and 
remains an independent variable in all the evolution equations. In order to use 
a physically meaningful variable in the planar case, one must convert all of 
the equations to a different variable, such as w = H -\- uj. 

Equation ((22|) allows us to assess how the planet appears to move from 
the point of view of a fixed observer: a planet will still appear to circulate 
around its parent star without ever reversing direction on the sky subject to 
all these external forces only if the planet is on a planar orbit. The extent of 
any modulation in the apparent motion, provided entirely by dfi/dt, increases 
as the planet's orbit becomes increasingly non-coplanar. 



4 Adiabatic Limit 

A planet on an orbit which is compact enough so that its period is smaller 
than the timescale of any external perturbations can be treated in the adia- 
batic approximation. This approximation has been utilized in many facets of 
planetary astronomy because the majority of exoplanets discovered are within 
hundreds of AU of their parent stars. 



4.1 Orbital Element Equations 

For the remainder of this paper, in order to simplify and focus our results, 
we assume Tuu — = '^ww = 0. These terms do not play a role in man y 
applications, such as tidal distortions of planetary orbits ( Brasser et al. . 2010l) . 



but were included in the general equations for completeness and future studies. 

Here, we consider the "adiabatic" case, where averaging over the fast vari- 
ables is justified. This corresponds to the following inequalities being satisfied: 
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T^xlr? < 1, T^yln^ < 1, Ty^^lr? < 1, Tyyjr?- < 1, T,,,Jr? < 1, T^^jn < 1, 
Tyu/n <C 1. We refer to this approximation as adiabatic, and this condition 
renders the equations of motion independent of the true anomaly /. We assume 



3/2 



dt 



df n (1 + e cos fY 



(23) 



which is the unperturbed two-body term. We obtain the adiabatic hmit by 
computing: 

dp _ v_ f'^'' dp 

dt adiabatic 27r 



-df 

dt non— adiabatic df 



(24) 



for a general variable /?. Carrying out the averaging, we find: 



da — 



cost 



dt n 
a 

~ 8? 



(-8 + + 8\/l - cos {2^2) sin (2w) cos ? 
2 + + 2a/1 - (3 + cos 2i) cos 2co 



- 2e^sin2«j sin 212 
de 5e\/l — 



(25) 



cos w sin LO sin zT^^ + 



dt 2n ' 16n 

+ [4 cos z cos 2lij sin 2J7 + sin 2a; (cos 2i7 (3 + cos 2i) + 2 sin^ z)] T^ra; 

— [4 cos i cos 2a; sin 2J7 + sin 2a; (cos 2J7 (3 + cos 2i) — 2 sin^ i)] Tyy 

— [4 cos z cos 2a; cos 2i7 — 4 cos i — sin 2a; sin 2 J? (3 + cos 2z)] T^y 



[4 cos % cos 2a; cos 2J7 + 4 cos i — sin 2a; sin 2i7 (3 + cos 2i)\ Tj 



yx 



+ 



8e3 



(1 - e^) (2 - - 2v/r^) 



X (4 cos i cos 2J7 sin 2a; + (3 + cos 2i) cos 2a; sin 2i7) 
dz 5e^ sin 2a; sin 2i 



{Tuv + J (26) 



-T 



smi 



dt %ns/Y^ 8nvT^ 

+ [sin 2f2 (2 + 3e^ + 5e^ cos 2a;) - lOe^ cos i sin 2a; sin^ D\ 
- [sin 2f? (2 + 3e^ + 5e^ cos 2a;) + lOe^ cos i sin 2a; cos^ /?] Tj 



+ [2 sin^ i? (2 + 3e^ + 5e^ cos 2a;) + lOe^ cos i sin 2a; cos fl sin 

— [2 cos^ i? (2 + 3e^ + 5e^ cos 2a;) — lOe^ cos i sin 2a; cos CI sin ]?] T^^: | 



+ - sin 2i sin 2i? (T„^ + T„„) 



(27) 
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~dt 



T 



cos z (2 + 3e^ - 5e^ cos 2a;) 

[siii^ J7 cos i (— 2 — 3e^ + 5e^ cos 2a;) + 5e^ cos J7 sin ]7 sin 2a;] 
[cos^ i7 cos ? (— 2 — 3e^ + 5e^ cos 2a;) — 5e^ cos f2 sin sin 2a;] Tj 
[cos J7 sin f2 cos z (— 2 — 3e 

[cos J7 sin n cos z (— 2 — 36^^ + 5e^ cos 2a;) + 5e^ cos^ i7 sin 2a;] Ty^ j- 



i sin^ fn Tuv + \ \ cos^ Q \ T„ 



(28) 



do; 



(sin^z-e^) - (1-e^) 
2nVl - e2 



T 



(2Cio + Cll + C12) 3^a:x + ( — 2Cio — Cii + C12) 

(— 2Cio cot 2J7 — Ci3 — lOe^ cos i sin 2a;) T^y 



2Cio cot 2J? — Ci3 + lOe^ cos z sin 2a;) T„; 



1 



(— 4Ci4 + C15) (T„^ + T^„) 



(29) 



16e4 

where the variables C'lo . . . C15 are given in Appendix A. When only vertical 
exter nal forces are included, the equations yield the same stationary solution 
as in iBrasseii (|2001[ ). In the more general case, when the non- vertical ext ernal 



forces are included, the equatio ns reduce to Eqs . (3)-(7 ) of FouchardI ( 2004h and 
the equations in Appendix A of Fouchard et al.l ( 20061 ). given their assumptions 
about T and after converting their Delaunay elements to the orbital elements 
used in this work. 

The Tzz, and Tyy terms in the expression for da/dt vanish. These terms 
are non-zero in all other orbital element equations. Therefore, a planet's semi- 
major axis is never secularly affected by vertical perturbations. In no case 
does orbit-averaging cause a velocity-dependent T term to vanish. However, 
in all cases, except for df2/dt, both velocity dependent terms appear in the 
combination (T„„ -t- T„„), which does vanish when there is no net shear on the 
two-body system. 



5 Planar Adiabatic Limit 

In the planar adiabatic limit, the general equations are greatly simplified. If 
the perturbations allow the planet to remain in its original orbital plane, then 
the vertical contributions vanish. The planet may have a planar orbit in two 
ways: in a prograde (z = 0°) or retrograde sense (z = 180°). As mentioned 
previously, because H is not defined physically for z = 0° or z = 180°, we wish 
to eliminate it from the equations of motion in the planar caseH. Trigonometric 

^ Mathematically, in the planar case, 17 7^ 0° and changes with time, as can be seen in 
Eq. {nil. 
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manipulation demonstrates that in both the prograde and retrograde cases, 
both n and lu will vanish when a suitable angle is defined. 



5.1 Prograde Equations 

In the prograde case, we use the longitude of pericenter vj = Q -\-uj and i = 0°. 
When the planet's orbit is sufficiently small, then we can impose the adiabatic 
approximation to obtain 



da 




'dt ^ 


n 


de 


5eVl - e2 


dt ~ 


4n 




(1 - e^) (2 



a ( - 2 + 2%/l - e^) sin 2vj 



U/IL. ^I^V-^ ^ ) ^^^^ ^ ^ 



[sin 2tu {T^x -T^yy) + (2 sin^ tu) T^y - (2 cos^ w) Ty^] 

- p2 _ 2^/l - e2) sin2n7 

^ (T„„+r„„) (31) 



2e3 



— — — [(3 + 5 cos2ri7) T-c^ + (3 — 5 cos 2ri7) Tyy + 5 s'm2Tu [T^y + '^yx)] 



-6e2-4(l-e2)"'"jcos2tn 



4e4 



(33) 



If {Tuv -\-T^u) 7^ 0, then the effect on the orbit may be drastic, even for orbiting 
"Hot Jupiters" at a 0.05 AU. The terms which represent the coefficients of 
(^uu + '^vu) for da/dt, de/dt and duj/dt are well-defined in the limits of e — > 
and e — )■ 1. Further, the term for de/dt takes on a maximum of w 0.06 at 

~ 0.68. This term is independent of a and n and hence might 
dominate the eccentricity evolution. 

Note too that the pericenter will be constantly perturbed by the external 
force regardless of any properties of the planet unless T^v = T^u = 0; even if 
Tuv = 7^ and e = 0, the last term in Eq. p3p will immediately become 
non-zero. 



5.2 Retrograde Equations 

The equations of motion for the completely retrograde case (i — 180°) may 
take on a similar form if a different angle, an obverse of pericenter, is defined 
to be the dog-leg angle equal to the sum of the longitude of descending node 
and the angle between the radius vector of the descending node and the peri- 
center of the orbit (i; = 180° — J7 -1-180° +a; — oj — n). Defining this angle allows 
us to eliminate both i? and uj, as in the previous subsection. Both nodes may 
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be physically important; the longstanding prevalence of the ascending node 
in dynamical parlance is perhaps due to the abundance of prograde orbits in 
the Solar System and, likely, exoplanetary systems, where prograde is defined 
with respect to the parent star's rotation. 

Thus, for the planar retrograde case, we find 



a\/l — e 



da 
'dt 



de 5eVT— ? 



[T^y - Ty^) + ^ + T,„) (34) 



di 



di 
dt 



An 



[sin2<j (T^x -^yy) - (2 sin^ <?) Ta;y + (2 cos^ <;) Yy^] 



(1 - e2) (2 - - 2%/!"^) sin2<; 



2e3 







(35) 
(36) 



d<; 



dt 



An 



[(3 + 5 cos 2(;) T^x + (3 - 5 cos 2<j) Tyy - 5 sin 2? {Y^y + Ty^) 



+ 



1 



4 + e'* 



66^ 



4 1 



,3/2 



COS 2^ 



4e4 



(37) 



Note that there is at least one sign change in a term from the prograde case 
for each of the expressions for da/dt, de/ dt and d<;/dt. 



6 An Application: Galactic Tides 



One potential application of the equations is the modelling of Galactic tides. 
There has been previous work on the effect of the Galactic tide on t he Port 



Clou d , as this tide strongly affects the orbital evolution of cloud comets (jHeisler fc: Tremaine 



19861: iMatese fc Whitmanl . Il989l ). At the Sun's location in the Galaxy, the 



dominant component of the tide is caused by the Galactic disk and it acts in 
a direction perpendicular to the disk. For the planets in the Solar system, the 
Galactic tide is not generally important. A rough rule-of-thumb is that the 
precession due to tides Ptido ~ ^cxt/^pb where Pcxt is the orbital period of the 
host star in the Galaxy and Ppi is the orbital period of the plant around the 
star. For Jupiter, this gives Ptide ~ 10^^ yr, well in excess of a Hubble time. 
The timescale for the effects of the Galactic tide on the planets in the Solar 
system to manifest themselves is very long. 

However, the growing consensus that exoplanets are ubiquitous throughout 
the Galaxy, together with the discovery of wide-or bit («1000-2500 AU) planets 
(e.g. Kuzuhara et al. . 2011 : Luhman et al. I l201ll) . motivates the study of the 
effects of Galactic tides much more generally than has been done hitherto. 

If the host star moving on a circular orbit in the Galaxy, t hen the equations 
for the effects of Galactic tides on the satellite of a are given bv lHeisler fc Tremaine 
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( 19861 ) and Brasser et al. ( 2010h . In this picture, the frame of reference is cen- 
tered on the star, orbiting with the star but not rotating. The frame is non- 
inertial. Assuming a completely flat Galactic rotation curve, the only non- 
vanishing planar terms are 



Txx = cos (2i7Gt) 
Tyx = i^G^ sin(2i7Gt) 



Txy = X2g' sin (2/2Gi) 
Tyy — ~Qg^ cos (2]7g^) 



(38) 



where Qq is the circular frequency of the star with respect to the Galac- 
tic center. The only non- vanishing vertical term is Tzz = —4:t:Gpq, where 
Pq is the local Galactic matter density (including both the luminous matter 
and any contribution from dark matterjfj. As characteristic of the Galacto- 
cent ric distances p r obed by microlensing - exemplifie d by OGLE 2007-BL G- 
050 ([Batista et al.l . l2009t ) and OGLE-2003-BLG-235 (iBennett et all . |2006|) - 
we take the location of the host star to be 3 kpc from the Galactic Center. We 
make the simple assumption that the Galactic rotation curve is flat with am- 
plitude 220 kms~^. Thus, for a Galactocentric distance R, we have f^G — 220 
kms^^/i?. Further, at i? = 3 kpc, we set pc — 0. 65Mf7^pc~^ , self-consistently 



2P 

generating the rotation curve (see e.g., Eq. 2.2 of lEvanslll993[ ). 

At a distance of 3 kpc, wide-orbit planets may no longer be bound to 
their parent star. One way to quantify th is boundary is to c o mput e the Hill 
radius of the star, Rh- Equation (20) of Jiang fc Tremain"3 ( 2010f) suggests 
that Rh = [Gm*/(2i7^)]i/3. Using fie = 220kms"V^, we find Rh = 7.3 x 
10'*AU(i?/kpc)2/3. Therefore, for i? = 3 kpc, the HiU radius is at least 2 x 10^ 
AU, suggesting that any orbiting objects within that distance will remain 
bound. 

We now apply these tides to the wide-orbit planet WD 0806-661B b, Sedna- 
like objects, and distant scattered disk objects. This last example provides a 
demonstration of when the adiabatic approximation breaks down, and hence 
when Eqs. must be used instead of Eqs. (P5|) -([ ^ . For every curve 

on every plot, we integrated the equations of motion in both orbital elements 
and Cartesian elements to check that the results are equivalent. 



6.1 Wide-orbit planet WD 0806-661B b 



Luhman et al.l ()2011l ) discovered the substellar companion to WD 0806-661B 



using direct imaging. This procedure yields a projected separation but no 
direct information about the orbital properties, including the eccentricity. 
Hence we assume the semimajor axis is equal to the projected separation 
(ao = a w 2500 AU), and provide representative values for other orbital prop- 
erties: Co = 0.5, vjQ = f2o = 0°. The dependence of the orbital evolution on 
the mass of the planet is negligible. 



^ T here is a sign error in the corresponding equation, the third Eq. 3, in iBrasser et al.l 
ll201Cri . 
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Widest-Orbit Exoplanet 



0.75 




2000 4000 6000 8000 10 000 



Time /My r 

Fig. 1 The eccentricity evolution of the exoplanet with the widest-known orbit, WD 0806- 
661B b, due to Galactic tides. Assumed initial parameters are ag = 2500 AU, eo = 0.5, cjq = 
0°,Qo = 0°. 



In Fig. m we plot the eccentricity evolution of WD 0806-661B b for 9 
different values of io over 10 Gyr, a characteristic main sequence lifetime. When 
iQ > 42°, the eccentricity changes by over 0.1. When io > 71°, the eccentricity 
changes by over 0.2. These variations illustrate that orbital signatures can 
significantly depend on the main sequence age of a wide-orbit planet. 



6.2 Sedna-like bodies 



Sedna represents the most distant member of the S olar System ever ob served. 



and currently resides about 90 AU from the Sun ([Brown et al.l . 120041 ). close 
to perihelion. A recent orbital fit in terms of orbital elements for Sedna has 
been computed by the Jet Propulsion Laboratorjlf] and shows that Oc « 544 
AU, Cc ~ 0.859, LOc w 310.9° and f2c « 144.42° where the subscript "c" stands 
for "current". Sedna is inclined with respect to the ecliptic by about 12°. The 
ability for Sedna to retain unperturbed by other Solar System objects until 
the end of the Sun's main sequence lifetime is uncertain; during the Sun's 
post-main se quence evolution, over 6 Gyr from now, Sedna will then evolve 
significantly (|Veras fc Wvattl . I2OI2I) . 

Because of this uncertainty, and because the purpose of this paper is to 
present the general perturbed equations of motion, we model the evolution of 
several Sedna-like objects. We fix oq = Oc and Co = Cc but vary iq, loq and 
J7oj and run the simulations for 10 Gyr to show the full extent of the orbital 
change over a typical main sequence lifetime. We present our results in Fig. [21 
Panels a, b and c correspond to {ljq ~ ujc, ■f^o — i^c), (^^o = 225°, Hq = 0°) and 



5 As of September 26, 2012, from [iittp://ssd.jpl .nasa.gov/sbdb. cgi?sstr=Sedna[ 
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Sedna-like Bodies 




2000 4000 6000 8000 10 000 
Time /My r 




Time/Myr 

Fig. 2 The eccentricity evolution of the Sedna-Uke objects, all with ao = 544 AU and 
eo = 0.859. Panels a, b and c correspond to (ujo = ujc, 1^0 = ^c), (i^o = 225°, Ho = 0°) 
and {loq = ■'^0 = 0°)i respectively. The modulation of the curves in panel c is due to planar 
tides. 
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Scattered Disk Bodies 




100 200 300 400 

Time/Myr 




100 200 300 400 

Time/Myr 

Fig. 3 The eccentricity evolution of initially nearly circular (eo = 0.05) scattered disk 
bodies (panel a: ao = 3 X lO** AU; panel b: ao = 5 X 10* AU; panel c: ao = 7 X 10** AU) in 
systems inclined at i = 60° with respect to the Galactic plane. The initial true anomalies of 
the objects are given in the legends in the upper panel, and loq = = 0° in all cases. In 
the adiabatic regime, all curves on a given plot would be equivalent. The plot demonstrates 
how adiabaticity is broken. 
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Planar Transition to Non-Adiabaticity 
Prograde Retrograde 
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100 
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Fig. 4 A comparison of prograde planar and retrograde planar transitions to non- 
adiabaticity. Here, i = 0° so vertical tides vanish. The left panels represent the prograde case 
with ■ujQ = 0° and the right panels represent the retrograde case with ? = 0°. The semimajor 
axes sampled here are equal to those in Fig. [3] (3 X 10* AU for the a panels, 5 X 10* AU 
for the b panels, and 7 X 10* AU for the c panels). The plot demonstrates how retrograde 
planets are generally more stable and more adiabatic than their prograde counterparts. 



Perturbed Two-Body Problem 



19 



{u!q — Hq — 0°), respectively. These sets of values were chosen to show differ- 
ent evolutionary behaviours. In panel a, the eccentricity decreases by several 
hundredths; in panel b, the eccentricity increases by several hundredths; in 
panel c, the eccentricity changes only by a few ten-thousandths. The curves in 
panel c also exhibit a visually discernible modulation due to the planar tides. 



6.3 The breaking of adiabaticity 

Here we demonstrate how adiabaticity is broken by considering outer scattered 
disk/inner Oort cloud-like objects with a = 3 x 10^ AU - 7 x 10^ AU. We 
integrate the general equations (Eqs. fTHOIH) in Fig. [31 as well as the planar 
prograde and retrograde limits in Fig. |4l 

Figure [3] illustrates in panel a that the transition region between adiabatic- 
ity and non-adiabaticity occurs just beyond about 3 x 10^ AU. In panel b, at 
flo = 5 X 10* AU, the sinusoidal form of the eccentricity evolution is heavily 
modulated, particularly after 200 Myr have elapsed. In panel c, at oq = 7 x 10* 
AU, all of the objects have become unstable before 400 Myr have elapsed. Also 
note that the initially small eccentricity of the orbiting bodies do not protect 
them from large eccentricity variations. 

Because planar tides are weaker than vertical tides, the adiabatic approxi- 
mation should remain robust to larger distances in the planar case than in the 
non-planar case (Fig. [3]). Fig. 2] demonstrates that the difference may be over 
a factor of 2 in semimajor axis the prograde case and over a factor of 3 in the 
retrograde case. Also, Fig. U] demonstrates that retrograde objects are more 
stable than prograde objects. This result is intuitively sensible, given that in 
the planar restricted circular three-body problem, the Hill sphere is largest for 
a tertiary orbiting t he secondary in the opposite s ense of the secondary's orbit 
about the primary ( Valtonen fc Karttunenl . [20061 pp. 131-133). Thus, objects 



on retrograde orbits are not likely to be as dynamically excited as objects on 
prograde orbits, or ones with nonzero inclinations. 



7 Conclusions 

This work provides an analytical characterization of the perturbed two-body 
problem. In principle, the perturbations T may represent any force; the po- 
sition and velocity-dependent perturbations we chose here are particularly 
well-suited for modelling the effect of Galactic phenomena on single-planet 
exosystems. 

The main results are the compact form (Eq. 2]) , from which perturbative 
equations of motion may be derived analytically and quickly using an algebraic 
software progranlfl, the general non-adiabatic non-planar equations of motion 
in terms of orbital elements for several types of perturbations (Egs. fTill^n)) . and 



^ A Mathematica notebook with the necessary precomputed matrices is available from 
the authors. 
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the explicitly-expressed specific cases of these equations (Eqs.j^H-UHl Eas-ISOl 
[551 and Eas. lMll57|) . We have discussed applications to wide-separation planets 
as well as distant trans-Neptunian objects in our Solar system. 

Further, there are numerous applications of these equations in planetary 
dynamics. For example, a number of studies have investigated the id ea of a 
Galactic habitable zone ( Gonzalez et al.l . l20Qll : lLineweaver et al. ■ 12004) . These 



are the regions in a galaxy in which a host star can harbor terrestrial planets 
that can support life. Others have argued that perhaps the entire Galaxy is 
(or h as been) suitable for life, and the idea of zones is too simplistic ( Prantzo^ 



20081 ). There are many processes influencing such a complex topic as habit- 



ability, including chemical, geophysical and biological factors. However, a pre- 
requisite for life is the existence of long-lived, stable planetary orbits around 
a host star. In particular, the planet must remain largely unaffected by per- 
turbations from Galactic tides, spiral arms, bars and rings, giant molecular 
clouds and dark matter substructures. The host star may lie on a nearly cir- 
cular orbit within the Galactic disk - much like the Sun - or it may lie on 
a more eccentric orbit belonging to the Galactic disk, spheroid or bulge. A 
substantial task for dynamicists is to map out the regions of our Galaxy in 
which host stars can provide planetary orbits that are stable against a rich 
variety of perturbations. The formalism that we have developed in this paper 
can play an important role in this ambitious quest. 
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A Some Auxiliary Variables 

Here, we collect some cumbersome variables referred to in the main body of the paper. The 
variables Ci, . . . Cg in the general case of Eqs. I I14I I- II20I I are given by 



Ci = e cos u) + cos {/ + u)) 
C2 = e sin u) + sin {/ + 

C3 = cos i sin Q sin ( f + u)) — cos fl cos (/ + oj) 

C4 = cos i cos Q sin ( f + u)) + sin fl cos (/ + oj) 

C5 = (3 -I- 4e cos / + cos 2/) sin u) + 2{e + cos /) cos u> sin / 

Ce = (3 -I- 4e cos / + cos 2/) cos u) — 2(e + cos /) sin u) sin / 

C-j = (3 + 2e cos / — cos 2/) cos u) + sin uj sin 2/ 

Cg = (3 — cos 2/) sin uj — 2{e + cos /) cos to sin / 

C9 = (3 + 2e cos / — cos 2/) sin u) — cos u> sin 2/ 



The variables Cio, . ■ ■ C15 in the adiabatic limit of Eqs. II25I I- II29II are given by 
Cio = 5 (e^ - 2) sin (2u)) sin {2J7) cos i 

Cii = cos2i7 (1 -66^-5 (-3 + 2e'^) cos2tj - 10 cos 2i sin^ oj) 

C12 = 11 - 66^ + cos2aj (5 - lOe^) + 10cos2isin^ w 

C13 = sin2i7 (-1 -f 6e^ + 5 (-3 + 2e^) cos 2w + 10cos2isin^ tj) 




